"""
Script to simulate the diode responce and determine the diode parameters
like dEdX, L_dead, etc...
"""

import readchn
import math
import numpy
import matplotlib.pyplot as plt
from isotopes import decays
import scipy


# Load the run we're going to use to fit diode performance
A = readchn.esc_run(run_path='130802E2.ZIP')
# We need an initial energy scale guess, we have a method that provides that!
A.initial_find_peaks()
A.smart_guess_peaks()
data = A.sum_spec


def diode_responce(energy,energy_scale,energy_offset,dEdX,dead_thickness,gauss_sigma,n=1e5):
   """
   This creats a simulated MCA of alpha events on a diode with a known energy scale and offset
   """
   # The MCA basically has 1024 channels, initially with 0 in all channels
   hist_array = numpy.zeros(1024)
   
   # Energy offset caused by the dead layer needs to be added to the energy scale
   # otherwise the peaks all end up lower than the found peaks in data
   dead_offset = dead_thickness * dEdX # energy
   
   for theta in numpy.random.random(n):
      if theta >= 0.5:
         #skip it, alpha went the wrong way
         continue
      phi = math.acos(2.*numpy.random.random() - 1.)

      #The length the track is in the silicone
      L = (energy+dead_offset)/dEdX

      # the dead layer wont make detectable energy
      L_dead = dead_thickness/math.sin(phi)

      #Energy remaining to deposit in active layer
      E_dep = (L-L_dead)*dEdX
      if E_dep < 1.:
         continue
   
      #Smear with a gaussian
      E_dep_smear = numpy.random.normal(E_dep,gauss_sigma,1)[0]
      #chns = 1024              # chanels in MCA
      #energy_scale = 80.       #chn/MeV
      #energy_offset = 0.1      #MeV @ chn=1
   
      # Put in the virtual MCA
      #print int(math.floor(E_dep_smear*energy_scale+energy_offset))
      hist_array[int(math.floor(E_dep_smear*energy_scale+energy_offset))] += 1
   print 'size=',hist_array.size
   return hist_array


def sim_spec(parms):
   a210, a214, a218, energy_scale, energy_offset, dEdX, thick_dead, gauss_sigma = \
        sim_parms[0], sim_parms[1],sim_parms[2],sim_parms[3],sim_parms[4],sim_parms[5],sim_parms[6],sim_parms[7]
   po210 = diode_responce(decays['Po210'][0]['energy'], energy_scale, energy_offset, dEdX, thick_dead, gauss_sigma)
   po214 = diode_responce(decays['Po214'][0]['energy'], energy_scale, energy_offset, dEdX, thick_dead, gauss_sigma)
   po218 = diode_responce(decays['Po218'][0]['energy'], energy_scale, energy_offset, dEdX, thick_dead, gauss_sigma)
   
   return a210*po210 + a214*po214 + a218*po218

def sim_nll(sim_parms):
   sim = sim_spec(sim_parms)
   return numpy.sum( sim - data * numpy.log(sim) )

def sim_res(sim_parms):
   sim = sim_spec(sim_parms)
   return numpy.sum( (sim - data)**2 )


#---ooo000OOO000ooo------ooo000OOO000ooo------ooo000OOO000ooo------ooo000OOO000ooo---
#               Thats all our functions we need, here the script
#---ooo000OOO000ooo------ooo000OOO000ooo------ooo000OOO000ooo------ooo000OOO000ooo---

#energy_scale=A.energy_parms[0], energy_offset=A.energy_parms[1],dEdX=0.11, thick_dead=0.01, gauss_sigma=0.09
sim_parms = [.3, .4, .35, A.energy_parms[0], A.energy_parms[1], 0.09, 2., 0.1]


parm_bounds = [(0.,1.),(0.,1.),(0.,1.),(A.energy_parms[0]*1.1,A.energy_parms[0]*.9), (A.energy_parms[1]*0.9,A.energy_parms[1]*1.1),\
                (0.0001,1.), (0.0001,1.), (0.001,.5)]

opt_parms = scipy.optimize.minimize(sim_res, sim_parms, method='L-BFGS-B', bounds= parm_bounds, tol=0.005, options={'disp':True})
print opt_parms

plt.plot(A.chn_energy,A.sum_spec,'b')
plt.plot(A.chn_energy,sim_spec(sim_parms),'r')
plt.show()

